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Abstract. The CuO-planes of high-T c superconductors were found to consist of geometric stripes with 
alternating superconducting and antiferromagnetic areas. Here we will investigate the repulsive Hubbard 
model of striped clusters as a possible microscopic description of the superconducting elements. The focus 
of our attention lies on the superconducting properties. We report in agreement with the square Hubbard 
model a signal in the d :E 2_ y 2-channel and investigate its dependence on system size, cluster shape and 
interaction strength. 
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1 Introduction 

Short after the discovery of the high-T c superconductors JlJ , 
the Hubbard model was introduced as a generic descrip- 
tion of the CuO-planes on a microscopic level Q . Accord- 
ing to the Van Hove scenario we use an extension, the 
tt'-Hubbard model, to shift the Van Hove singularity in 
the density of states close to the Fermi energy The 
experimental result of striped domains |4j,g|] in the super- 
conducting CuO-planes inside the high-T c materials has 
inspired this study of striped Hubbard clusters. 

In order to understand superconductivity in the high- 
Tc cuprates on a macroscopic level, the high-T c glass model 
was introduced in 1987 HQ]- 

It was demonstrated, that the high-T c glass model in- 
cluding the tt'-Hubbard model as a microscopic descrip- 
tion of the striped superconducting domains is able to ex- 
plain several properties of the high-T c cuprates || , e.sr the 
d-wave symmetry of the superconducting phase fe.llO] or 
the pseudogap above T c in the density of states lllQl^l. 
Furthermore this combined high-T c glass and tt'-Hubbard 
model picture gives an intuitive description of the experi- 
mental puzzle, that different samples of the same material 
and same doping exhibit a nearly constant superconduct- 
ing transition temperature T c , yet the critical current den- 
sities vary from sample to sample |HJ . 

Hubbard clusters were already investigated with nu- 
merical algorithms for a number of different geometries 
and dimensions, e.g. the one-dimensional chains and lad- 
ders f|l4lil|Jl|,|l^,@, two-dimensional (21^ squares |l9|, 



and layered square systems [^,|25|,|26[]. The 
stripe instability was found theoretically within Hartree- 
Fock calculations applied to an extended Hubbard model 
|^7| , and was confirmed by a number of subsequent inves- 



tigations |28| . But up to now it is not known whether the 
pure Hubbard model exhibits striping, e.g., in the form 
of a phase separation. In the closely related 2D t — J 
model the occurrence of stripes is discussed controversially 
Here we study striped clusters of the Hub- 
bard model directly and do not examine the occurrence of 
stripes per se, but only the existence of superconductivity 
in striped Hubbard clusters. 

In a single striped domain we consider the tt'-Hubbard 
model, which is described in real space by 



W = TLkt 



n 



poi 



with the kinetic 



'Hkin = 
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and the potential part 



U 



pot 
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of the Hamiltonian. We denote the creation operator for 
an electron with spin a at site i with c\ a , the correspond- 
ing annihilation operator with c$ a , and the number oper- 
ator at site i with m >a = c\ a c i a . The hopping tij is only 
nonzero for nearest neighbors i,j (ijj = t) and next near- 
est neighbors (Uj — t'). Finally U is the interaction. Usu- 
ally we choose t' < to shift the Van Hove singularity in 
the density of states close to the Fermi energy for less than 
half filled systems ((n) < 1, where (n) = (n e ,t + n e,j.)/2 
and n ei0 - is the number of electrons with spin a). Through- 
out this paper we set n e ,f = n e „[ = n e and the energy unit 
as t — 1. Additionally we apply periodic boundary condi- 
tions both in x- and y-direction. 



2 Hubbard model and superconducting 
correlations 

Following [[35],|3(| we use the (vertex) correlation function 
(instead of the largest eigen value of the reduced two par- 
ticle density matrix) as an indicator of superconductivity, 
making use of the standard concept of off-diagonal long 
range order (ODLRO) |7). 

We concentrate here on the da;2_ y 2-wave symmetry 
(abbreviated as d-wave). The full two-particle correlation 
function is defined for the the d-wave symmetry by 



CrfO) = 7 £ 9S9S' (4, ]4+s^+r+&' A^+ta) 

i,5,5< 



(4) 



The vertex correlation function C\ (r) is the two-particle 
correlation function Cd{r) without the contributions of the 
single-particle correlations of the same symmetry [Q . For 
the d-wave the result is 

CY(r) = C d (r)-±J2g s g s ,Ci(i,r)C L (i + 6,i+r + 5') . 

i,S,6' 

(5) 

In equation (y) C a (i, r) = (cj a c i+r CT ) is the single-particle 
correlation function for spin a. The phase factors are gs, 
gs' = ±1 to model the d-wave symmetry, the number of 
lattice points is L and the sum 5 (resp. 6') is over all 
nearest neighbors. 

We averaged the vertex correlation function Cj(r) 
only in the large range regime of r, i. e. for the distances 
|r| > \r c \: 

cY^t E c d( r ) (6) 

r,\r\>\r c \ 

with the number L c of lattice points with |r| > \r c \. The 
qualitative behavior of our results (concerning the vertex 
correlation functions) is not influenced by our choice of 
\r c \ as long as we suppress the short range correlations 
(i.e. r c > 1.9). 

Evidence for ODLRO in the d x 2_ y 2 channel was al- 
ready found for the case of the square 2D tt '-Hubbard 



model [El, 39 



We report here the existence of ODLRO 
in the striped Hubbard model and investigate the influence 
of shape and interaction strength on the superconducting 
signal. 

Figure [l] shows the d-wave correlation functions (eq. 
(H) and (H)) as a function of the distance |r| between the 
pair creation and pair annihilation operators of both the 
vertex and the full correlation function for a striped sys- 
tem. Similar to results for square systems, we obtain huge 
correlations for the short range part and finite, positive 
values for (r) for large distances in the system (inlay 
of figure 0). Here it becomes obvious that the vertex cor- 
relation function is non-negative in the d-wave case. 

For a comparison we plot in figure |l| also the correlation 
function for the extended s-wave (xs-) symmetry for the 
nearest neighbors. This symmetry obeys the same formu- 
las as equations (§|) and only the phase factors gs and 
gs> both are set equal to 1. In contrast to the d-wave case, 



the xs-wave symmetry does not exhibit this long range 
behavior (inlay of figure [l]) . This is also in agreement with 
simulations for the square Hubbard model |2l| , ^o| . 



3 The PQMC-Method 

We calculate the ground state properties of the Hubbard 
model using the projector quantum Monte Carlo (PQMC) 
method [0^] . In this algorithm the ground state is pro- 
jected with 



1 



(7) 



from a test state ^Pt) with the projection parameter 8 and 
the normalization factor N j43|l . In order to perform this 
projection it is necessary to transform the many-particle 
problem into a single-particle problem. This is done in 
two steps, first the exponential of the Hamiltonian TL is 
decomposed into two separate parts, Tikir, and 7i po t, us- 
ing a Trotter Suzuki transformation [[L|,|42(] and second 
the interaction term is treated with a discrete Hubbard 
Stratonovich (HS) transformation, which leads to an ef- 
fective single-particle problem with additional fluctuating 
HS fields g|. 

We use the second order Trotter Suzuki transforma- 
tion, which reads as 



(8) 

where m is the number of Trotter slices and r = — . Here 

m 

a systematic error of order C(r 2 ) enters the calculations 
for finite m. 

The two parameters m and 8 influence the correct pro- 
jection of the ground state l^o) from the test wave function 
|%) in the PQMC algorithm 

In figure |^ we investigate the dependence of the ground 
state energy per site Eq/L and of the vertex correlation 
function C d ' P on the Trotter parameter m. Both, Eq/L 
and C d ' , level off for large m, indicating the convergence 
of the PQMC method. The results resemble similar PQMC 
simulations for the square 2D-tt '-Hubbard model [^8| , |4s| . 
1^.0 and the APEX-oxygen model §!. Like in these 
cases the vertex correlation function is here more sensitive 
to m than the ground state energy per site Eq/L (figure ^) . 

Figure || shows the influence of the projection para- 
meter 9 on the same observables. The results are again 
in good agreement with similar simulations for the square 
Hubbard model 

Due to the more rapid convergence of the ground state 
energy compared to the vertex correlation function the 
relative changes of Eq/L of figures |^ (a) and || (a) are 
significantly different compared to their (b) counterparts 
showing the vertex correlation function. This is also ex- 
pressed by the very different scales of the corresponding 
y-axes. 

Figure § (b) shows additionally to a 12x4 system the 
results for a twice as large 24 x 4 system. Convergence 
occurs here at higher values for 9. namely 8 > 16. Due to 
the sign problem (inlay of figure g| (a)) we were not able to 



2 



perform simulations for 6 > 16 . Similar effects occur also 
for PQMC simulations of the square Hubbard model fjol . 

Quantum Monte Carlo simulations are often plagued 
with the sign problem (0|H^]. The average sign {sign) 
enters the calculation for the expectation value (A) of an 
observable A by 

{A) E^±^)May) {A +) (AT) (Q) 



E<7.<7' °"0 



{sign) 



Here a and a' are configurations of the HS field, w((T, a') 
is their weight, and A((T, a') is the expectation value of 
A for er and a' [Q. Now w(a, a') can have both positive 
and negative values, thus when used in a Monte Carlo 
algorithm for a transition propability one uses the right 
hand side of equation (^). (A + ) and (A~) denote the sep- 
arate averages of HS-configurations a and a 1 with positive 
resp. negative weights w(a, a'). Generally speaking QMC 
simulations are only meaningful for {sign) close to 1. 

The average sign {sign) is known to decrease for in- 
creasing system size L, interaction strength U and projec- 
tion parameter 9 (see the inlay of figure y resp. and [pli 
|52| ) . But a small average sign leads among others to large 
statistical errors in the Monte Carlo process and renders 
the simulation results meaningless. 

From the above analysis of the dependence of the ground 
state energy and the correlation functions on m and 9 we 
conclude, that the PQMC simulations are converged for 
U = 2 when 9 > 8 in smaller systems and 9 > 16 in larger 



systems. A ratio r 



of the projection parameter 



and the Trotter parameter was found to be sufficient for 
a correct decomposition. Due to the sign problem there 
is only a small range of the parameters, where 9 can be 
chosen sufficiently large, so that the investigation of the 
vertex correlation function and its long range behavior is 
possible. This is similar to the case of the square Hubbard 
model [H. 

For smaller systems we performed also some simula- 
tions of the Hubbard model using the stochastic diagonal- 
ization J53 , 54 , figure [| They compare also favorably with 
their PQMC counterparts, which is an additional indica- 
tion that the PQMC performs correctly. 



4 Superconductivity in stripes 

We now investigate the dependence of the superconduct- 
ing properties of the striped Hubbard model on system 
size, shape and interaction strength. 

The geometry has a quite significant effect on the mag- 
nitude of the correlation functions. For increasing width 
L y of the stripes, figure the average long range part of 

W,p ■ 



always show a higher superconducting signal than square 
systems of the same size L. 

This is rather surprising when one takes into account 
that in striped systems, on average, the distances \r\ be- 
tween pair creation and pair annihilation operators are 
larger than in square systems with the same number of 
sites L. In our view there are two effects which may in- 
crease the superconducting correlations. First the anisotropy 
of L x and L y which leads to more finite size shells. Finite 
size shells refer to the energy levels of the free Hubbard 
clusters (U — 0). It is known, that other ways of introduc- 
ing additional shells, e.g., anisotropic hopping t x ^ t y [55 
p6| , or additional hopping to next nearest neighbors [40 
increase the superconducting correlations in the repulsive 
Hubbard model. 

In our view it is a second effect, the squeezing of the 
system in one dimension, that gives rise to these increased 
superconducting correlations. 

In contrast to the width L y , the length L x of the stripes 
is relatively insensitive to the height of the plateau, fig- 
ure ^. 

Another way to strengthen the superconducting corre- 
lations is to increase the (repulsive) Hubbard interaction 
U, as shown in figure 0. Here we present both common 
methods for analy zing superconductivity: the full corre- 
lation function |57|,^5lp8l and the vertex function |j35U3c| , 
Eij . The dotted lines in figure @(a) and (b) indicate the 
values of the full (vertex) correlation function in the case 
of no interaction U = 0. Figure 0(a) shows that the full 
correlation function increases for higher interactions. The 
vertex correlation function is zero for the case of no in- 
teraction (U = 0) and increases also monotonous for in- 
creased interaction strength. Due to the sign problem we 
were not able to perform simulations for an interaction 
strength U > 2.5, for the system size and filling shown in 
figure 

Thus within the range of parameters accessible by the 
PQMC method, the superconducting correlations are in- 
creasing for an increasing repulsive interaction strength 
U. Both full and vertex correlation function show this 
behavior. We want to note, that the vertex correlation 
function is much more sensitive to variations of U due to 
the substraction of the background of the single-particle 
correlation functions (figure 0) . Thus the vertex correla- 
tion function is the more appropriate observable to ana- 
lyze superconductivity in small Hubbard clusters. These 
results are similar to observations made for the BCS re- 
duced Hubbard model [b4|. 



5 Effective interaction in striped Hubbard 
clusters 



the vertex correlation function C d ' is decreasing signif- Duc to thc failure of the ufmal finite gize gcali in the 

-yV. P - 



icantly. The ratio between C d ' for a rectangular 12x4 
and a square 12 x 12 system is almost 3. In figure § we 



show C d : for both, square systems and the rectangular 
12 x L y systems from figure U as a function of the system 
size L — L x x L y . Here the rectangular shaped systems 



square 2D Hubbard model we introduced |40| , |23 
effective model, the BCS-reduced Hubbard model, to com- 
pare the superconducting correlations for different system 
sizes L. This failure is mainly caused by thc underlying 
shell structure of the free (U = 0) system ^o|, Hl4oll . 



3 



The BCS-reduced Hubbard model exhibits the same 
corrections to scaling as the Hubbard model, and has a 
well chosen interaction term, that produces superconduc- 
tivity with d-wave symmetry. We calculate for this model 
the same correlation functions as for the Hubbard model. 
The effective interaction strength J e ff is then chosen to 
give the same values for the correlation functions as the 
Hubbard model (for details see J5(J ) . From this we get a di- 
rect estimate of the superconducting interaction strength. 

The calculation of an effective interaction for the three 
band Hubbard model was used to identify the pairing 
mechanism for d-wave superconductivity in this model. 
The evidence of d-wave pairing in this case is based on 
symmetry arguments and exact diagonalization results of 
small clusters 

In the momentum space the BCS-reduced Hubbard 
model is described by the Hamiltonian: 



H 



BCS 



■xjBCS 



H 



BCS 



(10) 



The kinetic part is again equation (|2|) , only transformed to 
momentum space, i. e. = Hkin, and the interaction 

is denned by (for d-wave interaction) 



H 



BCS 



p u fe,T c -fc,i £ 



-p,l c p,T 



(11) 



In equation (pi] ) we use the form factors 

f k = cos(k x ) - cos(fcy) (12) 

to model the d-wave symmetry in 2D (k = (k x , k y )). 

We calculate the ground state of this BCS-reduced 
Hubbard model with the exact and the stochastic diag- 
onalization [153 ,63, 5 4f| - 

In figure g and 3| we show the effective interaction J e f f 
corresponding to the correlation functions and systems 
shown in figure |] and [| Within the error bounds of the 
simulations we conclude that J e ff is nearly constant for 
various geometries of the system. It is not possible to cal- 
culate stochastic errors of the physical observables within 
the SD. But for smaller system sizes our comparison of 
SD with exact diagonalization results indicates that for 
the weak interactions J used here the errors in the SD are 
negligible J|0 , 54 . The error bars shown in figures f| and 
H are therefore calculated using only the statistical errors 
of the PQMC results and fitting these values to the SD 
results. 

In addition to the above mentioned, one has to take 
into account, that even so we tried to perform the calcula- 
tions at a constant filling (n) w 0.8 the constraint of closed 
shells for PQMC simulations leads to different fillings (n) 
for each of these system sizes. Furthermore, in the case of 
figure §| (and [| respectively) one has to take into account, 
that all simulations are performed at a constant 9 = 8. 
Whereas figure ^ indicates that for large system sizes L 
a higher value of 8 would lead to slightly higher values 
of the vertex correlation function in the PQMC runs and 
thus to a slightly lower effective interactions J e ff- 



From figures |8j and |9| we conclude that within the ac- 
curacy of the applied methods, the effective interaction 
strength J e // is equal for both square and striped sys- 
tems. Furthermore J e ff is insensitive to the length of the 
striped systems. 



6 Summary and Conclusions 

Here we performed ground state simulations of the 2D tt'- 
Hubbard model and the BCS-reduced Hubbard model of 
striped clusters using PQMC and SD techniques. Together 
with the exact diagonalization these are the most reliable 
computational tools for this type of calculations. 

We concentrated our investigations on the behavior of 
rectangular striped systems. In agreement with previous 
calculations of the square Hubbard model we find that 
these finite systems show evidence for superconductivity 
in the d x 2_ y 2 channel for repulsive interactions U. Com- 
pared to the squared case these correlations are signifi- 
cantly enhanced, and the superconducting signal is nearly 
insensitive to the length of these stripes. 

Using SD-techniques we were capable of estimating the 
effective superconducting interaction strength J e ff of a 
BCS-reduced Hubbard model with the same symmetry of 
the superconducting correlation functions. Within the ac- 
curacy of our methods both square and striped Hubbard 
model show approximately the same superconducting in- 
teraction strength J e ff- 

In conclusion, the striped Hubbard model is a promis- 
ing candidate for the microscopic description of the su- 
perconducting striped domains in the high-T c cuprates. 
Within the larger framework of the high-T c glass model a 
combined model is able to explain many puzzling proper- 
ties of the high-T c materials. 
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lull 
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Fig. 1. Distance |r| between pair creation and pair annihila- 
tion operators versus the two-particle resp. vertex correlation 
function of the Hubbard model with U = 2. System parame- 
ters: L x = 12, L a = 4, n e = 21, i' = -0.22, 6 = 8, m = 64. 
The points labeled with /uZ/ and vertex are correlation func- 
tions with d-wave symmetry and the points labeled with xs 
show the vertex correlation function with xs-symmetry. 
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Fig. 3. ^-scaling. System parameters: L x = 12, L y = 4, n e = 
15, J7 = 2, t' = -0.22, r = 1/8 (solid lines), and L x = 24, 
L„ = 4, n e = 41, U = 2, t' = -0.22, r = 1/8 (dashed line), 
(W = 1.9). 
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Fig. 4. Averaged vertex correlation function Cj' P for in- 
creasing L y . System parameters: L x = 12, U = 2, t' = —0.22, 
6» = 8, m = 64, and fillings: L y = 4: (n) = 0.88, = 6: 
(n) = 0.86, L H = 8: (n) = 0.77, L y = 10: (n) = 0.78 and 
L y = 12: (n) = 0.85. 



Fig. 2. m-scaling. System parameters: L x — 12, L y — 4, 
n e = 15, U = 2, t' = -0.22, = 8, full line |r c | = 1.9, dashed 
line |r c | = 2.9. 
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Fig. 5. Averaged vertex correlation function C d ' for in- 
creasing L y compared with square system sizes L — L x ■ L x . 
Simulation parameters: U = 2, t' = —0.22, 6 = 8, m = 64. 1.) 
Striped systems: (solid lines), sizes L x , L y , and fillings (n) see 
figure 4. 2.) Square systems: (dashed lines), sizes and fillings: 
L x = 6: (n) = 0.72, L x = 8: (n) = 0.78, = 10: (n) = 0.82 
and L x = 12: (n> = 0.85. 




Fig. 7. Averaged vertex (|r c | > 1.9) and full (all lattice points) 
d-wave correlation function for increasing interaction U. Sys- 
tem parameters: L = 12 x 4, (n) = 0.88, t' = -0.22, = 8, 
m — 64. The inlay shows the average sign. 
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Fig. 6. 

ing L x 



Averaged vertex correlation function C d ' P for increas- 
System parameters: L y = 4, U = 2, t' = -0.22, 9 = 8, 
m = 64. PQMC runs are averaged for r > r c = 1.9 (o) and 
l r l > l r c| — 2.9 (x). SD runs are averaged for \r\ > \r c \ = 1.9 
(©). Sizes and fillings: = 6: (n) = 0.92, L x = 8: (n) = 0.81, 
L x = 10: (n) = 0.75, L x = 12: (n) = 0.88, = 16: (n) = 0.78, 
L x = 20: (n) = 0.83 and L x = 24: (n) = 0.85. 



Fig. 8. Effective d-wave interaction J e /y for various L y . 
tern parameters: = 12, U = 2, i' = -0.22, 9 = 8, 
m = 64. Fillings see figure 4. 
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Fig. 9. Effective d-wave interaction J e ff for various L x . Sys- 
tem parameters: L y = 4, U = 2, t' = —0.22, 9 = 8, and 
m = 64. Fillings see figure 6. 
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